{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "16ade37b-0d5e-49bc-968e-57b36b60390f",
   "metadata": {},
   "outputs": [],
   "source": [
    "import cantera as ct\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy\n",
    "\n",
    "# parameter values\n",
    "p0 = 101325  # pressure\n",
    "T0 = 296  # inlet temperature\n",
    "phi = 0.7  # equivalence ratio\n",
    "alpha = 0.566  # volume fraction (<1)\n",
    "\n",
    "# mdot(=mass flux, 단위 면적당 시간당 전달되는 질량)\n",
    "mdot_reactants = 80  # kg/m^2/s # 반응물이 1초당 1m² 면적을 통해 80kg 만큼 들어온다\n",
    "mdot_products = mdot_reactants  # kg/m^2/s # 생성물이 반응 후에 같은 속도로 배출된다, mass conservation\n",
    "\n",
    "# Mechanism and fuel composition\n",
    "xh2 = alpha / (1 - alpha) * 1  # xh2 = alpha/(1-alpha)*xc\n",
    "fuel = {'CH4': 1, 'H2': xh2}\n",
    "oxidizer = {'O2': 1, 'N2': 3.76}\n",
    "\n",
    "# Calculation params\n",
    "width = 0.2  # m domain width\n",
    "loglevel = 0  # amount of diagnostic output (0 to 5)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "bc624fa9-4b5e-4360-9a92-5d0184a46b0f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Flame speed: 0.31479 m/s\n",
      "1812.8710116169764\n"
     ]
    }
   ],
   "source": [
    "# Gas object 1 for the counterflow flame\n",
    "gas1 = ct.Solution('gri30.yaml')  # reaction mechanism file\n",
    "gas1.TP = T0, p0\n",
    "\n",
    "# fix the composition of the fuel\n",
    "gas1.set_equivalence_ratio(phi, fuel, oxidizer)  # hold temperature and pressure constant\n",
    "\n",
    "# calculate flame speed using laminar unstretched flame\n",
    "flame = ct.FreeFlame(gas1, width=width)\n",
    "flame.solve(loglevel, auto=True)\n",
    "flame_speed = flame.velocity[0]  # flame velocity(first one)\n",
    "print(f'Flame speed: {flame_speed:.5f} m/s')\n",
    "\n",
    "# Create the counterflow premixed flame simulation object\n",
    "fl1 = ct.CounterflowPremixedFlame(gas=gas1, width=width)\n",
    "fl1.transport_model = 'multicomponent'\n",
    "fl1.energy_enabled = True  # energy equation\n",
    "fl1.set_refine_criteria(ratio=3, slope=0.1, curve=0.2, prune=0.02)\n",
    "\n",
    "# set the boundary flow rates\n",
    "fl1.reactants.mdot = mdot_reactants\n",
    "fl1.products.mdot = mdot_products\n",
    "\n",
    "fl1.set_initial_guess()  # assume adiabatic equilibrium products\n",
    "fl1.solve(loglevel, auto=True)\n",
    "print(fl1.products.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "20fe6486-fbbf-4893-9a5a-7618321036cf",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Index of max |du/dx|: 10\n",
      "Value of max |du/dx|: 2786.6408924401885\n"
     ]
    }
   ],
   "source": [
    "# calculate velocity gradient\n",
    "grad_u = np.gradient(flame.velocity, flame.grid) \n",
    "\n",
    "# find index that has max|du/dx|\n",
    "max_grad_idx = np.argmax(np.abs(grad_u))\n",
    "print(f\"Index of max |du/dx|: {max_grad_idx}\")\n",
    "\n",
    "# value of max|du/dx|\n",
    "max_grad_value = grad_u[max_grad_idx]\n",
    "print(f\"Value of max |du/dx|: {max_grad_value}\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "e258fd3d-2043-47d2-b6f0-5e7a12cf8536",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Thermal diffusivity of reactants (α): 2.56978e-05 m^2/s\n"
     ]
    }
   ],
   "source": [
    "# calculate thermal diffusivity\n",
    "cp = flame.cp_mass[0]  # specific heat capacity at constant pressure(Cp), first grid point (J/kg/K)\n",
    "density = flame.density_mass[0]  # density at the first grid point (kg/m^3)\n",
    "lambda_ = flame.thermal_conductivity[0]  # thermal conductivity at the first grid point (W/m/K)\n",
    "\n",
    "thermal_diffusivity = lambda_ / (density * cp)  # thermal diffusivity calculation\n",
    "print(f\"Thermal diffusivity of reactants (α): {thermal_diffusivity:.5e} m^2/s\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "0993a255-c13b-445f-8c92-a0e45415da8c",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Strain rate (Ka): 2.92952e-08 1/s\n"
     ]
    }
   ],
   "source": [
    "# calculate strain rate\n",
    "Ka = thermal_diffusivity / (flame_speed * max_grad_value)\n",
    "print(f\"Strain rate (Ka): {Ka:.5e} 1/s\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "1d771899-50ef-48fd-96f4-84b057809c2a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkkAAAHFCAYAAADmGm0KAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjkuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8hTgPZAAAACXBIWXMAAA9hAAAPYQGoP6dpAABzU0lEQVR4nO3dd1hT5/sG8DuBEIaAskEQcIsTRevGAbir1VattnXX1r1qtePraOuedde66q5Wra0LHLhARdyCuMCFiCICCkIg7+8PS35GooIGMrg/1+WlefPmnOfJgNtzTs6RCCEEiIiIiEiNVNcFEBEREekjhiQiIiIiDRiSiIiIiDRgSCIiIiLSgCGJiIiISAOGJCIiIiINGJKIiIiINGBIIiIiItKAIYmIiIhIA4YkKlY++ugjWFhY4MmTJ6+d07NnT8hkMjx48CDfy5VIJJg4ceL7F5jP5UdFRWHixImIi4vT6nomTpwIiUSi+mNmZgZvb28MHz78jc/Zu4iLi0O7du1gZ2cHiUSCESNGIC4uDhKJBKtXr1bNW716NSQSidZ7LQqLFy9W60VbmjVrBolEgrJly0LTRROOHDmieg0LY/1vM3HiRHh5eeV7rkQiea/15eTkwMnJCXPnzn3jPC8vL/Tu3fud1uHl5VWon3HSTwxJVKz069cPz58/x4YNGzTen5KSgu3bt6N9+/ZwdnYu4upeLzw8HP3791fdjoqKwqRJkwotOOzduxfh4eHYtWsXOnXqhAULFqBNmzYafyG/q5EjR+LkyZNYuXIlwsPDMXLkSLi6uiI8PBzt2rXT2np0qbBCEgBYW1sjNjYWBw8ezHPfypUrYWNjUyjr1UdHjhzBw4cP0blzZ12XQkbGVNcFEBWlNm3awM3NDStXrsSgQYPy3L9x40ZkZGSgX79+Oqju9erXr1+k66tTpw4cHBwAAIGBgUhKSsLatWsRFhaGRo0aaXxMeno6LC0t872OS5cuoV69eujUqZPaeFH3aqjKlCkDa2trrFy5Ei1btlSNp6WlYcuWLejZsyeWL1+uwwqLztatW+Hn5wdPT09dl0JGhluSqFgxMTFBr169EBkZiYsXL+a5f9WqVXB1dUWbNm0AAAkJCRg4cCDc3d1Vu54mTZqE7Ozst67r0qVL6NixI0qVKgVzc3PUqlULa9asyTPvyZMnGD16NMqWLQu5XA4nJye0bdsWV65cUc15eXfb6tWr8cknnwAAmjdvrrZb5aeffoKpqSnu3LmTZz19+/aFvb09nj9/nq/n6mW5weXWrVsAXuzuqVatGo4cOYKGDRvC0tISffv2BQDcvn0bn332GZycnCCXy1GlShXMnj0bSqUSABAaGgqJRILr169jz549qvrj4uI07m57nf3796Nly5awsbGBpaUlGjVqhAMHDrz1ca/bfZdbV2hoqGost8+jR4+ifv36sLCwQOnSpfHjjz8iJyfnjevx8vLC5cuXcfjwYVWPL++CetvzlB99+/bFtm3b1HaFbtq0CQDQvXv3PPOvX7+OPn36oEKFCrC0tETp0qXRoUOHPJ+Fr776Cubm5oiMjFSNKZVKtGzZEs7Ozrh//36+a3zZrl27UKtWLcjlcnh7e2PWrFl55rzpPaBpt7YQAtu3b0eXLl1UYwqFAmPHjoWLiwssLS3RuHFjnDp1Ks/j2rZtC3t7e9y+fVs1np6ejqpVq6JKlSp49uzZO/VJxoMhiYqdvn37QiKRYOXKlWrjUVFROHXqFHr16gUTExMkJCSgXr162LdvH/73v/9hz5496NevH6ZOnYoBAwa8cR0xMTFo2LAhLl++jF9//RXbtm2Dj48PevfujRkzZqjmpaWloXHjxli2bBn69OmDf/75B0uXLkXFihVf+4uoXbt2mDJlCgBg0aJFCA8PV+2iGjhwIExNTbFs2TK1xzx+/BibNm1Cv379YG5uXuDn7Pr16wAAR0dH1dj9+/fx2WefoUePHti9ezcGDRqEhw8fomHDhggODsZPP/2EnTt3IiAgAGPGjMGQIUMAALVr10Z4eDhcXFzQqFEjVf2urq75rmfdunUICgqCjY0N1qxZgz///BN2dnZo1apVvoJSQSQkJKB79+7o2bMn/v77b3z88cf4+eefMXz48Dc+bvv27Shbtix8fX1VPW7fvh0A8vU85Uf37t1hYmKCjRs3qsZWrFiBjz/+WOPutvj4eNjb22PatGnYu3cvFi1aBFNTU3zwwQeIiYlRzZs3bx6qVKmCrl27qgLYpEmTEBoainXr1hXotcp14MABdOzYEdbW1ti0aRNmzpyJP//8E6tWrSrwsl4WFhaG+/fvq4WkAQMGYNasWfjiiy/w999/o0uXLujcuTOSk5NVcyQSCdauXQtLS0t07doVCoUCADBo0CDExsbizz//hJWV1XvVRkZAEBVD/v7+wsHBQWRlZanGRo8eLQCIq1evCiGEGDhwoChRooS4deuW2mNnzZolAIjLly+rxgCICRMmqG53795dyOVycfv2bbXHtmnTRlhaWoonT54IIYSYPHmyACBCQkLeWO+ry9+yZYsAIA4dOpRnbq9evYSTk5PIzMxUjU2fPl1IpVIRGxv7xvVMmDBBABAJCQlCoVCI5ORksW7dOmFhYSE8PDxERkaGEOLF8wdAHDhwQO3x48aNEwDEyZMn1ca//vprIZFIRExMjGrM09NTtGvXTm1ebGysACBWrVqlGlu1apUAoKr92bNnws7OTnTo0EHtsTk5OaJmzZqiXr16b+zx1eXlOnToUJ7nNLfPv//+W23ugAEDhFQqzfPeeFXVqlWFv79/nvGCPE+a+Pv7i6pVqwohXrzefn5+QgghLl++LACI0NBQERERkee5fFV2drbIysoSFSpUECNHjlS779q1a8LGxkZ06tRJ7N+/X0ilUvHDDz+8sa5cEyZMEJ6enmpjH3zwgXBzc1O9h4QQIjU1VdjZ2YmXfxVpeg/kevVzIIQQI0aMENWrV1fdjo6OFgDy9LN+/XoBQPTq1Utt/NixY8LU1FSMGDFCrFy5UgAQv//+e551e3p65lk3GT9uSaJiqV+/fnj06BF27twJAMjOzsa6devQpEkTVKhQAQDw77//onnz5nBzc0N2drbqT+6uuMOHD792+QcPHkTLli3h4eGhNt67d2+kp6cjPDwcALBnzx5UrFgRAQEBWutt+PDhSExMxJYtWwC82E2yZMkStGvXLt/fOHJxcYFMJkOpUqXw2WefoXbt2ti7d6/aVqhSpUqhRYsWao87ePAgfHx8UK9ePbXx3r17Qwih8SDjggoLC8Pjx4/Rq1cvtddFqVSidevWiIiI0OpuEmtra3z44YdqYz169IBSqcSRI0feaZnafJ769u2L06dP4+LFi1ixYgXKlSuHpk2bapybnZ2NKVOmwMfHB2ZmZjA1NYWZmRmuXbuG6Ohotbnly5fH8uXLsWPHDrRv3x5NmjR55293PXv2DBEREejcubPae8ja2hodOnR4p2Xm2rZtm9pWpEOHDgF48S3Vl3Xt2hWmpnkPw23UqBF++eUXzJs3D19//TU+++wzvTsmkXSHIYmKpY8//hi2traqTf27d+/GgwcP1H44PnjwAP/88w9kMpnan6pVqwIAHj169NrlJyUladwl4ebmprofeLHbxd3dXWt9AYCvry+aNGmCRYsWAXgR9uLi4gq0G2f//v2IiIjAuXPn8OjRIxw7dgw+Pj5qczT1l9++30fuqRk+/vjjPK/N9OnTIYTA48eP33s9uTR9y9HFxQXAu/ejzeepadOmqFChApYtW4a1a9eqdidrMmrUKPz444/o1KkT/vnnH5w8eRIRERGoWbMmMjIy8sxv164dnJ2d8fz5c4waNQomJib5rutlycnJUCqVquftZZrG8uvUqVO4ffu2WkjKfe5eXa6pqSns7e01Lqdnz54wMzNDZmYmvvnmm3euh4wPv91GxZKFhQU+/fRTLF++HPfv38fKlSthbW2tOiAaABwcHFCjRg388ssvGpeR+wtNE3t7e43HFMXHx6uWDbw4xufu3bvv04pGw4YNwyeffIIzZ85g4cKFqFixIgIDA/P9+Jo1a6pqfB1Nv4jz2/f7yF3GggULXvtNuDedviF3S0ZmZqba+OtCr6bzZSUkJADAa3/pvo22n6c+ffrghx9+gEQiQa9evV47b926dfjiiy9Ux7TlevToEUqWLJln/ldffYW0tDRUrVoVw4YNQ5MmTVCqVKkC1Qa82OookUhUz9vLXh173eujKTj+9ddfqFixIqpVq6Yay31NEhISULp0adV4dna2xmXk5OSgZ8+eKFWqFORyOfr164fjx4/DzMysAB2SseKWJCq2+vXrh5ycHMycORO7d+9G9+7d1b7C3r59e1y6dAnlypWDn59fnj9vCkktW7bEwYMHVb/0cv3xxx+wtLRU/XJv06YNrl69WuDdUHK5HAA0/u8feHHSzDJlymD06NHYv38/Bg0a9N4n7MuPli1bIioqCmfOnFEb/+OPPyCRSNC8efP3XkejRo1QsmRJREVFaXxd/Pz83vgLLneX44ULF9TGc3e9viotLS3PfRs2bIBUKn3tbq1ccrlc42uk7eepV69e6NChA7755hu1YPAqiUSieu/k2rVrF+7du5dn7u+//45169Zh4cKF2LlzJ548eYI+ffoUqK5cVlZWqFevHrZt26b27cq0tDT8888/anOdnZ1hbm6e5/X5+++/8yz3r7/+UtuKBLz4RiIArF+/Xm38zz//1Pit1AkTJuDo0aNYv349Nm/ejPPnz3NrEv0/HR8TRaRTNWrUEBKJRAAQJ06cULsvPj5eeHp6isqVK4vFixeLAwcOiF27dolFixaJdu3aiTt37qjm4pUDSq9cuSKsra1FxYoVxbp168Tu3btFz549BQAxY8YM1bzU1FRRtWpVUaJECfHzzz+L4OBg8ffff4tRo0aJgwcPvnb5N2/eFABEp06dxNGjR0VERIR49OiRWv3Tp08XAISVlZXqQPG3yT1w++HDh2+c9/KBwy9LTEwUpUuXFi4uLuK3334T+/btE8OGDRMSiUQMGjRIbe67HrgthBBr164VUqlUdOvWTWzZskUcPnxYbN26Vfz444/iq6++emPt2dnZolKlSqJMmTJiw4YNYs+ePeLLL78U3t7eGg/ctre3F25ubmLBggVi3759Yvjw4QKA+Prrr9+4HiFeHFQtl8vFpk2bxKlTp8SFCxcK/Dxp8rrn/2WaDtz+4osvhFwuF3PnzhUHDhwQM2bMEI6OjsLd3V3tAPMLFy4ICwsLtYOct27dKgCIuXPnvrU+TQduBwcHC6lUKho3biy2b98utm7dKurWrSs8PDzEq7+K+vfvL8zNzcXs2bPF/v37xZQpU0S1atXUPgdnz54VAMTp06fzrP+zzz4TEolEjB07VgQHB4s5c+YINzc3YWNjo9ZTbk0vf7Zyv5ixbds2tWXywO3iiSGJirX58+cLAMLHx0fj/Q8fPhTDhg0T3t7eQiaTCTs7O1GnTh3x/fffi6dPn6rmvRpihBDi4sWLokOHDsLW1laYmZmJmjVravzGTnJyshg+fLgoU6aMkMlkwsnJSbRr105cuXLljcufN2+e8Pb2FiYmJhq/DRQXFycAvDU0vOx9Q5IQQty6dUv06NFD2NvbC5lMJipVqiRmzpwpcnJy1Oa9T0gSQojDhw+Ldu3aCTs7OyGTyUTp0qVFu3btxJYtW97a59WrV0VQUJCwsbERjo6OYujQoWLXrl0aQ1LVqlVFaGio8PPzE3K5XLi6uorvvvtOKBSKt64nLi5OBAUFCWtrawFALTjk93nS5F1DUnJysujXr59wcnISlpaWonHjxuLo0aPC399fFZKePn0qKleuLHx8fMSzZ8/Uljl48GAhk8nyfCvvVZpCkhBC7Ny5U9SoUUOYmZmJMmXKiGnTpqnecy9LSUkR/fv3F87OzsLKykp06NBB9X7O/Rz88MMPGtchhBCZmZli9OjRwsnJSZibm4v69euL8PBw4enpqQpJ8fHxwsnJSbRo0ULtOVcqlaJDhw6iZMmSau85hqTiSSKEFq8zQER6Y8GCBRg2bBguXbqkOticCqZZs2Z49OgRLl26pOtSDMrEiROxevXqQr3eno+PD9q0aYPZs2cX2jpelnvdN16/rXjhgdtERubs2bOIjY3F5MmT0bFjRwYkMkpRUVG6LoGKAYYkIiPz0UcfISEhAU2aNMHSpUt1XQ4RkcFiSCIyMoW5i6O4efk6bkRU/PCYJCIiIiINeJ4kIiIiIg0YkoiIiIg04DFJ70CpVCI+Ph7W1tZFchZjIiIien9CCKSlpcHNzQ1S6du3EzEkvYP4+Pg8V3cnIiIiw3Dnzp18XVycIekdWFtbA3jxJNvY2GhtuQqFAsHBwQgKCoJMJtPacvWFsfcHGH+Pxt4fYPw9sj/DZ+w9FmZ/qamp8PDwUP0efxuGpHeQu4vNxsZG6yHJ0tISNjY2RvvGN+b+AOPv0dj7A4y/R/Zn+Iy9x6LoL7+HyvDAbSIiIiINGJKIiIiINGBIIiIiItKAxyQREZFW5eTkQKFQ6GTdCoUCpqameP78OXJycnRSQ2Ez9h7fpz+ZTAYTExOt1cKQREREWiGEQEJCAp48eaLTGlxcXHDnzh2jPY+dsff4vv2VLFkSLi4uWnluGJKIiEgrcgOSk5MTLC0tdfILXKlU4unTpyhRokS+ThZoiIy9x3ftTwiB9PR0JCYmAgBcXV3fuxaGJCIiem85OTmqgGRvb6+zOpRKJbKysmBubm6UAQIw/h7fpz8LCwsAQGJiIpycnN5715vxPbtERFTkco9BsrS01HElVNzlvge1cVwcQxIREWmNMR4jQ4ZFm+9BhiQiIiIiDRiSiIiI3pOXlxfmzZunt8szFKGhoTAxMUFKSspr56xevRolS5YsknoYkoiIqNjq0KEDAgICNN4XHh4OiUSCM2fOFHFVQEREBL788kvVbYlEgh07dhR5He8iISEBw4cPR/ny5WFubg5nZ2c0btwYS5cuRXp6+hsf27BhQ9y7d0+r10V9H/x2mx65lZSOlCwgPSsbNqam3LdPRFTI+vXrh86dO+PWrVvw9PRUu2/lypWoVasWateuXeR1OTo6Fvk6teHmzZto1KgRSpYsiSlTpqB69erIzs7G1atXsXLlSri5ueHDDz/U+FiFQgEzMzO4uLggNTW1iCvXjFuS9EjHxeH4X6Qpav50EOW/34Oak4LRaNpBtJ53BB8vCUOfVacwdONZjN92EVN3R2PBgWtYdTwWWyPvYu+lBBy//ggX7j5BQspzCCF03Q4Rkd5r3749nJycsHr1arXx9PR0bN68Gf369QMAhIWFoWnTprCwsICHhweGDRuGZ8+evXa5t2/fRseOHVGiRAnY2Niga9euePDggdqcnTt3ws/PD+bm5nBwcEDnzp1V9728u83LywsA8NFHH0EikaBs2bK4ffs2TE1Ncfr0abVlLliwAJ6enhp/B4wfPx7169fPM16jRg1MmDABwIvdXfXq1YOVlRVKliyJRo0a4datW6/t81WDBg1S1dW1a1dUqVIF1atXR5cuXbBr1y506NBBNVcikWDp0qXo2LEjrKys8PPPP2vc3bZ69WqUKVMGlpaW+Oijj5CUlJTvet4XtyTpCaVSQGYihQTZEJAgRymQkqFASsa7fYXR2twUFZ2tUdG5BCo4Wav+7Wgt5xYqIioSQghkKIr2shlKpRIZWTmwzud/FE1NTfHFF19g9erV+N///qf6+bhlyxZkZWWhZ8+euHjxIlq1aoWffvoJK1aswMOHDzFkyBAMGTIEq1atyrNMIQQ6deoEKysrHD58GNnZ2Rg0aBC6deuG0NBQAMCuXbvQuXNnfP/991i7di2ysrKwa9cujTVGRETAyckJq1atQuvWrSGRSCCXy9GyZUusWrUKfn5+qrmrVq1C7969Nf6c79mzJ6ZNm4YbN26gXLlyAIDLly/j4sWL2Lp1K7Kzs9GpUycMGDAAGzduRFZWFk6dOpXv3xlJSUkIDg7GlClTYGVlpXHOq8uaMGECpk6dirlz58LExASxsbFq9588eRJ9+/bFlClT0LlzZ+zdu1cV6IoCQ5KekEoliPiuOXbt2o3mgUHIzJEgLTMbac+z8fR5Np5mKpCq+nc20p4r/vs7+///fv5i/EFaJtKeZyPyVjIibyWrrcfWQoYKTiVQ4b/QVNHZGhWcS8CxBMMTEWlXhiIHPv/bp5N1X5oYiBL5PJFg3759MXPmTISGhqJ58+YAXuxq69y5M0qVKoXhw4ejR48eGDFiBACgQoUK+PXXX+Hv748lS5bA3NxcbXn79+/HhQsXEBsbCw8PDwDA2rVrUbVqVURERKBu3br45Zdf0L17d0yaNEn1uJo1a2qsL3fXW+7lNpRKJVJTU9GvXz8MGjQIc+bMgVwux/nz53Hu3Dls27ZN43KqVauGGjVqYMOGDfjxxx8BAOvXr0fdunVRsWJFPH78GCkpKWjfvr0qRFWpUiVfzyEAXL9+HUIIVKpUSW3cwcEBz58/BwAMHjwY06dPV93Xo0cP9O3bV3X71ZA0f/58tGrVCuPGjQMAVKxYEWFhYdi7d2++63ofDEl6RiIBLM1MYSuTwekdl5GZnYPYR89w9cFTXHuQhqsP0nDtwVPEJT1DSoYCp28l4/Qr4amkpQwVnV4EporO1qjpURLVS9vCRMrgRETGrXLlymjYsCFWrlyJ5s2b48aNGzh69CiCg4MBAJGRkbh+/TrWr1+veowQAkqlErGxsXmCRHR0NDw8PFQBCQB8fHxQsmRJREdHo27dujh37hwGDBjwXnV36tQJw4YNw/bt29G9e3dV/bm75zTp2bMnVq5ciR9//BFCCGzcuFEV/uzs7NC7d2+0atUKgYGBCAgIQNeuXQt8eY9X/8N96tQpKJVK9OzZE5mZmWr3vbwVTJPo6Gh89NFHamMNGjRgSKJ3Jzc1QWUXG1R2Uf92wHNFDm4+fIZrif8fnK4lPsWtpGd4kq7AqbjHOBX3WDW/pKUMTSo4omkFB/hXdISTjfmrqyIiei0LmQmiJrcq0nUqlUqkpabBQlawy1H069cPQ4YMwaJFi7Bq1Sp4enqiZcuWqmUOHDgQw4YNy/O4MmXK5BkTQmjcMv/yeO7lM96HmZkZPv/8c6xatQqdO3fGhg0b3nragB49emDcuHE4c+YMMjIycOfOHXTv3l11/6pVqzBs2DDs3bsXmzdvxg8//ICQkBCNxzK9qnz58pBIJLhy5YraeNmyZQFo7vl1u+Vy6fr4WoakYsRcZgIfNxv4uOUNTzcePsW1B09x9UEaYhLScCruMZ6kK/DP+Xj8cz4eAFDZxRr+lRzhX8ERdbxKQW76ftfEISLjJpFIYGlWtL9mlEolss1MCnz4QNeuXTF8+HBs2LABa9aswYABA1TLqF27Ni5fvozy5cvna1k+Pj64ffs27ty5o9qaFBUVhZSUFNVWpxo1auDAgQPo06dPvpYpk8mQk5P3+K7+/fujWrVqWLx4MRQKhdrB35q4u7ujadOmWL9+PTIyMhAQEABnZ2e1Ob6+vvD19cX48ePRoEEDbNiwIV8hyd7eHoGBgVi4cCGGDh361gCUHz4+Pjhx4oTa2Ku3CxNDEsFcZoKqbrao6marGsvOUeLcnSc4fPUhjlx9iAv3UnAlIQ1XEtKw7PBNWJqZoEFZe/hXckTTCo7wcnj/DwMRka6UKFEC3bp1w3fffYeUlBT07t1bdd+3336L+vXrY/DgwRgwYACsrKwQHR2NkJAQLFiwIM+yAgICUKNGDfTs2RPz5s1THbjt7++v2r00YcIEtGzZEuXKlUP37t2RnZ2NPXv2YOzYsRrr8/LywoEDB9CoUSPIZDLVhVurVKmC+vXr49tvv0Xfvn3ztYWqZ8+emDhxIrKysjB37lzVeGxsLH777Td8+OGHcHNzQ0xMDK5evYovvvgCwIvdZl988QUOHDiA0qVLa1z24sWL0ahRI/j5+WHixImoUaMGpFIpIiIicOXKFdSpU+et9b1s2LBhaNiwIWbMmIFOnTohODi4yHa1ATwFAL2GqYkUfl52GB1UCX8PaYzIHwIxv3stdK5dGg4l5EjPysGBK4n439+X0WxWKPxnHsKPOy5hf9QDPMvM1nX5REQF1q9fPyQnJyMgIEBtN1qNGjVw+PBhXLt2DU2aNIGvry9+/PHH1x6rk3vix1KlSqFp06YICAhA2bJlsXnzZtWcZs2aYcuWLdi5cydq1aqFFi1a4OTJk6+tbfbs2QgJCYGHh0eeoNGvXz9kZWWpHQD9Jp988gmSkpKQnp6OTp06qcYtLS1x5coVdOnSBRUrVsSXX36JIUOGYODAgQBenBYhJibmjReOLVeuHM6ePYuAgACMHz8eNWvWhJ+fHxYsWIAxY8bgp59+yleNuerXr4/ff/8dCxYsQK1atRAcHIwffvihQMt4HxKh6x1+Big1NRW2trZISUnR6llBFQoFdu/ejbZt20Imk2ltudqmVApEJ6TiyNVHOHw1EZG3kqHI+f+3kcxEAj9PO7Sp7oKONUvD1vJFL4bS3/sw9h6NvT/A+HssrP6eP3+O2NhYeHt75/m2V1HK/eaXjY0NpFLj3A7wao+//PILNm3ahIsXL+q6NK1439fwTe/Fgv7+5u42KjCpVKLaPfd1s3J4mpmN8BtJOHw1EYevPsSdxxkIv5mE8JtJ+GVXNNpWd0W3uh6o7W6t69KJiIzG06dPERMTgwULFhR4Cw3lj8HG7KlTp0Iikai+ugi8OAp+4sSJcHNzg4WFBZo1a4bLly+rPS4zMxNDhw6Fg4MDrKys8OGHH+Lu3btFXL1xKSE3RaCPM37uVB1HvmmOQ2Oa4Yd2VVDZxRqZ2UpsP3sP3X87gaD5x7H/ngQP0zLfvlAiInqjoUOHonHjxvD398/3rjYqGIMMSREREfjtt99Qo0YNtfEZM2Zgzpw5WLhwISIiIuDi4oLAwECkpaWp5owYMQLbt2/Hpk2bcOzYMTx9+hTt27fX+K0BKjiJRAJvByv0b1IWe4Y3wd+DG+HTemVgZWaCuKR0/HPbBE1nHcHAtadx6EoicpTc20tE9C5WrVqFzMxMbN68WXUgN2mXwYWkp0+fomfPnli+fDlKlSqlGhdCYN68efj+++/RuXNnVKtWDWvWrEF6ejo2bNgAAEhJScGKFSswe/ZsBAQEwNfXF+vWrcPFixexf/9+XbVktCQSCWp6lMTUztVx6vsATOlUFV4lBLKVAvsuP0Cf1RFoNO0g5gTH4M7jN18ZmoiIqKgZ3DFJgwcPRrt27RAQEICff/5ZNR4bG4uEhAQEBQWpxuRyOfz9/REWFoaBAwciMjISCoVCbY6bmxuqVauGsLAwtGql+aRnmZmZamcJzb06sUKheONR/gWVuyxtLlNfmEmBTjWcYPUgB141G2LHhQfYcS4eCanP8evB61hw6DoalrVH1zql0bKKE+SmBpffARj3awgYf3+A8fdYWP1lZ2dDCIGcnBwolUqtLrsgcr+LlHtGbGNk7D2+b385OTkQQiA7OzvP+7yg73uDCkmbNm3CmTNnEBERkee+hIQEAMhzUixnZ2fVFYwTEhJgZmamtgUqd07u4zWZOnWq2vV1cgUHB8PS0rLAfbxNSEiI1pepT+LOh6EWgGrVgYuPJQhPlCAmRYrjN5Jw/EYSrEwF6joKNHBSwkX7T2+RMPbX0Nj7A4y/R233J5FI4OrqisePH8PaWvdf0nj5MAtjZew9vmt/aWlpePbsGQ4ePJjnjN3p6QXba2EwIenOnTsYPnw4goOD3/j10lfPsvq608MXZM748eMxatQo1e3U1FR4eHggKChI66cACAkJQWBgoNF+9fjV/j787747yenYGhmPv87cw4O0TITelyD0vhQNy9lheIvyqF2mpM7qLoji+BoaG2PvsTD7e/DgAVJTU2Fubg5LS0udXDRbCIFnz57BysrKaC/abew9vmt/Qgikp6cjLS0Nrq6uqFWrVp45uXuC8stgQlJkZCQSExPVTqKVk5ODI0eOYOHChYiJiQHwYmvRyyf4SkxMVG1dcnFxQVZWFpKTk9W2JiUmJqJhw4avXbdcLodcLs8zLpPJCuWHaGEtV19o6q+sky3GtrHFqKBKOHLtITaduoMDVxIRduMxwm6cQpMKDhgRUBF1PEu9Zqn6pTi+hsbG2HssjP5Kly4NExMTPHr0SKvLLQghBDIyMmBhYWGUAQIw/h7ft79SpUrBxcVF42ML+p43mJDUsmXLPCfK6tOnDypXroxvv/0WZcuWhYuLC0JCQuDr6wsAyMrKwuHDhzF9+nQAQJ06dSCTyRASEoKuXbsCAO7fv49Lly5hxowZRdsQaWRqIkWLys5oUdkZdx6nY9Gh69gaeRdHrz3C0WuP0LSiI0YGVIBvGcMIS0TFSe4uNycnJ50d06VQKHDkyBE0bdrUaEOusff4Pv29fMkWbTCYkGRtbY1q1aqpjVlZWcHe3l41PmLECEyZMgUVKlRAhQoVMGXKFFhaWqJHjx4AAFtbW/Tr1w+jR4+Gvb097OzsMGbMGFSvXh0BAQFF3hO9mYedJaZ1qYHBzctj4cHr2HrmLo78dy25ZpUcMSKgImp5lNR1mUT0ChMTE519Jd3ExATZ2dkwNzc3ygABGH+P+tSfwYSk/Bg7diwyMjIwaNAgJCcn44MPPkBwcLDaQYRz586FqakpunbtioyMDLRs2RKrV6/mOSb0mIedJaZ//CIsLTh4DdvO3kNozEOExjxEi8pOGBFQATXcS+q6TCIiMjIGHZJCQ0PVbkskEkycOBETJ0587WPMzc2xYMECjVduJv1Wxt4SMz+p+WLL0qHr2H72Hg5eScTBK4loWdkJIwIqorq7ra7LJCIiI2GYJ6OhYs3LwQqzPqmJA6P80bl2aUglwIErieiw8Bj6rzmNS/dSdF0iEREZAYYkMlheDlaY07UW9o/yx0e+L8LS/ugHaL/gGAb8cRqX4xmWiIjo3TEkkcEr61gCc7vVQsgof3Sq5QaJBAiJeoB2vx7DwLWncT3RuE+4RkREhYMhiYxGOccSmNfdFyEj/fFhzRdhad/lB2gz/yhm7ruCjCxexJiIiPKPIYmMTnmnEvj1U1+EjGyKlpWdoMgRWHToBgLnHsbBKw90XR4RERkIhiQyWuWdrPF7Lz8s+7wO3GzNcTc5A31Xn8ZXayMR/yRD1+UREZGeY0gioyaRSNCqqgtCRvnjy6ZlYSKVYO/lBATMOYzlR25CkWN8V9AmIiLtYEiiYsFKborv2lbBrmGN4edZCulZOfhldzQ6LDiGyFuPdV0eERHpIYYkKlYqu9jgz4ENML1LdZS0lOFKQhq6LAnHuL8uIPlZlq7LIyIiPcKQRMWOVCpBt7plcHB0M3xSxx0AsCniDlrOOYwtp+9ACKHjComISB8wJFGxZWdlhpmf1MSWrxqgonMJPH6WhW+2XkC3ZSdw9QHPrUREVNwxJFGxV9fLDruGNcG4NpVhITPBqbjHaDv/KKbtuYL0rGxdl0dERDrCkEQEQGYixVf+5RAyqikCfZyRrRRYevgGAuccQUgUz61ERFQcMSQRvcS9lCWWf+GH5V/4oXRJC9x7koEBf5zG8E1nkfZcoevyiIioCDEkEWkQ6OOMkFFNMdD/xbmV/j4Xj3a/HsO5O090XRoRERURhiSi17A0M8X4NlXw58D6KF3SArcfp+PjJWFYevgGlEp+A46IyNgxJBG9RR1PO+we3gTtqrsiWykwbc8VfLHyFBJTn+u6NCIiKkQMSUT5YGshw8IevpjepTrMZVIcu/4IbeYfxaEriboujYiICglDElE+SSQvTkL579DGqOxijaRnWeizOgI//RuFzOwcXZdHRERaxpBEVEDlnayxY3Aj9G7oBQBYcSwWnReH4ebDp7otjIiItIohiegdmMtMMPHDqvj9Cz+UspThcnwq2i84hr/O3AOvakJEZBwYkojeQ4CPM/YMb4oGZe2RnpWDcdsv449rUp5TiYjICDAkEb0nF1tzrOv/Ab5pVQkmUgnOJEnx4eITOHs7WdelERHRe2BIItICE6kEg5uXx8Z+dWEnF7ibnIFPloZjceh1nlOJiMhAMSQRaZFvmZIYWyMH7aq5IFspMGNvDD5feZLnVCIiMkAMSURaZmEKzO1aHTO61ICFzATHryeh9fyjOH79ka5LIyKiAmBIIioEEokEXet64J+hjeHjaoPHz7LwxcpTWH/ylq5LIyKifGJIIipE5Z1KYNughuhUyw05SoHvt1/CpH8uIztHqevSiIjoLRiSiAqZucwEc7vVwpigigCAVcfj0G/NaaTyNAFERHqNIYmoCEgkEgxpUQFLetaGuUyKw1cfosviMNxOStd1aURE9BoMSURFqE11V2wZ2BDONnJcS3yKTouP41TsY12XRUREGjAkERWx6u622DmkMaqXtsXjZ1no+fsJbDl9R9dlERHRKxiSiHTA2cYcfw5sgLbVXaDIEfhm6wVM3RONHJ54kohIbzAkEemIhZkJFn5aG8NalAcALDt8EwPXRuJZZraOKyMiIoAhiUinpFIJRgVVwvzutWBmKsX+6Af4eGk47j3J0HVpRETFHkMSkR7oWKs0Ng6oD4cSZoi+n4qOC4/jDC+QS0SkUwxJRHqijmcp7BjcCJVdrPHoaSa6/3YCf5+7p+uyiIiKLYYkIj3iXsoSW79uiIAqTsjKVmL4pnOYE3IVSh7QTURU5BiSiPRMCbkpln3uh4FNywIAfj1wDUM3nUVGVo6OKyMiKl4Ykoj0kIlUgvFtq2DGxzUgM5Fg14X76PZbOBLTnuu6NCKiYoMhiUiPdfXzwLp+H6CUpQwX7qag69Jw3HnMS5kQERUFhiQiPfdBWXtsH9QI7qUsEJeUjk+WhuPagzRdl0VEZPQYkogMgJeDFbZ+1RAVnEogIfU5ui4Lx/k7T3RdFhGRUWNIIjIQLrYvLmVS090WyekK9Fh+AmE3Hum6LCIio8WQRGRASlmZYf2A+mhYzh7PsnLQe1UEgi8n6LosIiKjxJBEZGBKyE2xsnddBPk4Iytbia/Xn8G2M3d1XRYRkdFhSCIyQOYyEyzuWRtdarsjRykw6s/zWH08VtdlEREZFYYkIgNlaiLFzI9roE8jLwDAxH+iMH//NQjBs3MTEWkDQxKRAZNKJfhfex+MCqwIAJi7/yom/xvFy5gQEWkBQxKRgZNIJBjWsgImdvABAKw6Hodvtl5Ado5Sx5URERk2hiQiI9G7kTfmdK0JE6kEf525i2GbzkLBoERE9M4YkoiMSOfa7lj6WR2YmUix+2ICRmw6x6BERPSOGJKIjEygjzOWfl77xYVxL97HiM3nuOuNiOgdMCQRGaEWlZ2x9LM6L4LSBQYlIqJ3wZBEZKRaVnHGkp4vgtK/F+5j5J/nGZSIiAqAIYnIiAX4OGNRjxe73v45H4/RWxiUiIjyiyGJyMgFVXXBwh61YSqV4O9z8Riz5TxyeB4lIqK3YkgiKgZavRSUdjAoERHlC0MSUTHRupoLFvbwhYlUgu1n7+EbBiUiojdiSCIqRlpXc8XCT18EpW1n72Hs1gsMSkREr8GQRFTMtKnuil+7+6rOzP3tXxd4rTciIg0YkoiKoXY1XDG/ey2YSCXYGnkX47ddhBAMSkREL2NIIiqm2tdww7xuL4LS5tN3MGV3NIMSEdFLGJKIirEONd0wvUsNAMDyo7FYHHpDxxUREekPhiSiYu7jOu74oV0VAMDMfTFYf/KWjisiItIPDElEhP5NymJI8/IAgB92XMK/F+J1XBERke4ZTEiaOnUq6tatC2trazg5OaFTp06IiYlRmyOEwMSJE+Hm5gYLCws0a9YMly9fVpuTmZmJoUOHwsHBAVZWVvjwww9x9+7domyFSC+NDqqInh+UgRDAyM3ncPjqQ12XRESkUwYTkg4fPozBgwfjxIkTCAkJQXZ2NoKCgvDs2TPVnBkzZmDOnDlYuHAhIiIi4OLigsDAQKSlpanmjBgxAtu3b8emTZtw7NgxPH36FO3bt0dOTo4u2iLSGxKJBJM7VkP7Gq5Q5Ah8tTYSkbeSdV0WEZHOGExI2rt3L3r37o2qVauiZs2aWLVqFW7fvo3IyEgAL7YizZs3D99//z06d+6MatWqYc2aNUhPT8eGDRsAACkpKVixYgVmz56NgIAA+Pr6Yt26dbh48SL279+vy/aI9IKJVII5XWuhaUVHZChy0Hd1BGIS0t7+QCIiI2QwIelVKSkpAAA7OzsAQGxsLBISEhAUFKSaI5fL4e/vj7CwMABAZGQkFAqF2hw3NzdUq1ZNNYeouDMzlWLpZ7VRu0xJpGQo8PmKk7idlK7rsoiIipyprgt4F0IIjBo1Co0bN0a1atUAAAkJCQAAZ2dntbnOzs64deuWao6ZmRlKlSqVZ07u4zXJzMxEZmam6nZqaioAQKFQQKFQvH9D/8ldljaXqU+MvT/AeHqUSYBlPX3Rc0UEriY+xWcrTmBT/3ooaf7i/1WG3t+bGMtr+Drsz/AZe4+F2V9Bl2mQIWnIkCG4cOECjh07luc+iUSidlsIkWfsVW+bM3XqVEyaNCnPeHBwMCwtLfNZdf6FhIRofZn6xNj7A4ynx888gPkpJrj9OAMfLwzF0Ko5sDQ1nv7exNh7ZH+Gz9h7LIz+0tMLtlXc4ELS0KFDsXPnThw5cgTu7u6qcRcXFwAvtha5urqqxhMTE1Vbl1xcXJCVlYXk5GS1rUmJiYlo2LDha9c5fvx4jBo1SnU7NTUVHh4eCAoKgo2NjdZ6UygUCAkJQWBgIGQymdaWqy+MvT/AOHts2DQdny4/hfinWdiSYI9PXZPQrrXx9PcqY3wNX8b+DJ+x91iY/eXuCcovgwlJQggMHToU27dvR2hoKLy9vdXu9/b2houLC0JCQuDr6wsAyMrKwuHDhzF9+nQAQJ06dSCTyRASEoKuXbsCAO7fv49Lly5hxowZr123XC6HXC7PMy6TyQrlDVpYy9UXxt4fYFw9lne2xR/9PkC3ZeE4cycFz9OkaNfWxGj6ex1jeg01YX+Gz9h7LIz+Cro8gzlwe/DgwVi3bh02bNgAa2trJCQkICEhARkZGQBe7GYbMWIEpkyZgu3bt+PSpUvo3bs3LC0t0aNHDwCAra0t+vXrh9GjR+PAgQM4e/YsPvvsM1SvXh0BAQG6bI9Ir1VxtcHK3nVhLpMi6okUE//hdd6IyPgZzJakJUuWAACaNWumNr5q1Sr07t0bADB27FhkZGRg0KBBSE5OxgcffIDg4GBYW1ur5s+dOxempqbo2rUrMjIy0LJlS6xevRomJiZF1QqRQfLzssP8bjXx1boz+DPyHjwdSmDwf2fpJiIyRgYTkvLzv1aJRIKJEydi4sSJr51jbm6OBQsWYMGCBVqsjqh4aFHJEV28ldgaa4KZ+2LgXsoCHWuV1nVZRESFwmB2txGRfmjiItC3oScA4JstFxAR91jHFRERFQ6GJCIqsG9bVUSrqs7IylFiwB+ncfPhU12XRESkdQxJRFRgUqkE87r5oqZHSTxJV6DP6ggkPc18+wOJiAwIQxIRvRMLMxP8/oUfPOwscCspHQP+OI3nCl4omoiMB0MSEb0zR2s5VvWuCxtzU5y5/QSj/zwPpZKnBiAi48CQRETvpbyTNZZ97geZiQS7Lt7HjH0xui6JiEgrGJKI6L01KGePGR/XAAAsPXwD60/e0nFFRETvjyGJiLTiI193jAyoCAD439+XERqTqOOKiIjeD0MSEWnNsJbl0aW2O3KUAoPXn0FUfMEuJklEpE8YkohIayQSCaZ2ro6G5ezxLCsHfVdH4H5Khq7LIiJ6JwxJRKRVZqZSLPmsDio4lUBC6nMMXBvJUwMQkUFiSCIirbO1kGFl77ooZSnDhbspGL/tYr6uv0hEpE8YkoioUHjYWWJRz9owkUqw/ew9/H40VtclEREVCEMSERWahuUc8L/2PgCAqXuicfjqQx1XRESUfwxJRFSovmjgiW5+HlAKYOiGM4h99EzXJRER5QtDEhEVKolEgsmdqqKOZymkPs9G/zURSHuu0HVZRERvxZBERIVObmqCJZ/VhouNOW48fIYRm87xGm9EpPcYkoioSDhZm+O3L+pAbirFgSuJmBNyVdclERG9EUMSERWZGu4lMa1LdQDAwkPX8e+FeB1XRET0egxJRFSkPvJ1x5dNywIAvtlyAZfjU3RcERGRZgxJRFTkvm1dGU0qOCBDkYMv/4hE0tNMXZdERJQHQxIRFTkTqQQLP60Nbwcr3HuSga/Xn4EiR6nrsoiI1DAkEZFO2FrKsPyLOighN8Wp2MeY/E+UrksiIlLDkEREOlPeyRrzutWCRAKsPXELG0/d1nVJREQqDElEpFMBPs4YE1QJADDh78u4cPeJbgsiIvoPQxIR6dygZuUQ6OOMrBwlvl53Bk/Ss3RdEhERQxIR6Z5EIsGsT2qijJ0l7j3JwOg/z/OM3ESkcwxJRKQXbC1kWNyzNsz+OyP3ksM3dF0SERVzDElEpDeqlbbF5A+rAgBmB8cg7MYjHVdERMUZQxIR6ZVudT3wcR13KAUwbONZPEh9ruuSiKiYYkgiIr0ikUjwU8dqqOxijUdPszB0w1meaJKIdIIhiYj0joWZCRb3rP3iRJNxjzFrX4yuSyKiYoghiYj0UlnHEpjxcQ0AwLIjN7HvcoKOKyKi4uadQtKdO3dw9OhR7Nu3D2fOnEFmJi9OSUTa17a6K/o19gYAjNlyHreSnum4IiIqTvIdkm7duoXx48fDy8sLXl5e8Pf3R5s2beDn5wdbW1sEBgZiy5YtUCp57AARac+4NpVRx7MU0p5n4+t1Z/BckaPrkoiomMhXSBo+fDiqV6+Oa9euYfLkybh8+TJSUlKQlZWFhIQE7N69G40bN8aPP/6IGjVqICIiorDrJqJiQmYixcIevrCzMkPU/VRM3HlZ1yURUTFhmp9JZmZmuHHjBhwdHfPc5+TkhBYtWqBFixaYMGECdu/ejVu3bqFu3bpaL5aIiidXWwvM714LX6w8hU0Rd1DHsxQ+8fPQdVlEZOTyFZJmzpyZ7wW2bdv2nYshInqdJhUcMTKgIuaEXMUPOy6hqpstfNxsdF0WERkxfruNiAzGkObl0bSiIzKzlRi0PhJPM7N1XRIRGbF8bUl6mbe3NyQSyWvvv3nz5nsVRET0OlKpBPO61UL7X48iLikdP2y/iLndar3xZxIR0bsqcEgaMWKE2m2FQoGzZ89i7969+Oabb7RVFxGRRnZWZvj1U190++0EdpyLR6PyDjw+iYgKRYFD0vDhwzWOL1q0CKdPn37vgoiI3sbPyw4jAypgVvBV/O/vy/AtUwrlnUrouiwiMjJaOyapTZs2+Ouvv7S1OCKiN/q6WXk0LGePDEUOhmzg+ZOISPu0FpK2bt0KOzs7bS2OiOiNTP47PsneygxXEtIwZXe0rksiIiNT4N1tvr6+agdJCiGQkJCAhw8fYvHixVotjojoTZxszDG7a030XhWBP8JvoWE5B7Su5qLrsojISBQ4JHXq1EnttlQqhaOjI5o1a4bKlStrqy4ionxpVskJA/3LYtnhmxi79TyqlbaBeylLXZdFREagwCFpwoQJhVEHEdE7GxNUCSdvPsa5O08wYtM5bPqyPkxNeBo4Ino/+fop8uxZwa68XdD5RETvQ2YixYJPfWEtN8XpW8lYdOiGrksiIiOQr5BUvnx5TJkyBfHx8a+dI4RASEgI2rRpg19//VVrBRIR5YeHnSV+6lQNAPDrwWuIvJWs44qIyNDla3dbaGgofvjhB0yaNAm1atWCn58f3NzcYG5ujuTkZERFRSE8PBwymQzjx4/Hl19+Wdh1ExHl0cm3NEJjErHjXDyGbzqL3cObwMZcpuuyiMhA5SskVapUCVu2bMHdu3exZcsWHDlyBGFhYcjIyICDgwN8fX2xfPlytG3bFlIpjwMgIt2Z3KkaIm8n487jDPxvxyXM6+6r65KIyEAV6MBtd3d3jBw5EiNHjiyseoiI3ouNuQzzuvmi67Jw7DgXj2aVnNDJt7SuyyIiA8TNPkRkdOp4lsKwFhUAAD/suIQ7j9N1XBERGSKGJCIySoObl4OfZyk8zczG8E1nkZ2j1HVJRGRgGJKIyCiZmkgxt1stWMtNceb2E54WgIgKjCGJiIyWh50lfv7o/08LcPY2TwtARPnHkERERq1jrdL4sKYbcpQCIzefw7PMbF2XREQGosAhycvLC5MnT8bt27cLox4iIq37qWM1uNqaIy4pHT/vitJ1OURkIAockkaPHo2///4bZcuWRWBgIDZt2oTMzMzCqI2ISCtsLWWY3bUmJBJg46k72B/1QNclEZEBKHBIGjp0KCIjIxEZGQkfHx8MGzYMrq6uGDJkCM6cOVMYNRIRvbeG5RzQv7E3AGDctgt49JT/uSOiN3vnY5Jq1qyJ+fPn4969e5gwYQJ+//131K1bFzVr1sTKlSshhNBmnURE7210UCVUcrbGo6dZGL/tIn9OEdEbvXNIUigU+PPPP/Hhhx9i9OjR8PPzw++//46uXbvi+++/R8+ePbVZJxHRezOXmWBut1qQmUgQEvUAWyPv6rokItJjBbosCQCcOXMGq1atwsaNG2FiYoLPP/8cc+fOReXKlVVzgoKC0LRpU60WSkSkDT5uNhgZWBEz9sZg0j9RaFDOHu6lLHVdFhHpoQJvSapbty6uXbuGJUuW4O7du5g1a5ZaQAIAHx8fdO/eXWtFEhFp08Cm/3827m+2XIBSyd1uRJRXgUPSzZs3sXfvXnzyySeQyWQa51hZWWHVqlXvXRwRUWEwkUowu2tNWJqZIPxmElaHxem6JCLSQwUOSc2bN0dSUlKe8SdPnqBs2bJaKYqIqLB52lvhu7ZVAADT917BjYdPdVwREembAoekuLg45OTk5BnPzMzEvXv3tFIUEVFR6PlBGTSp4IDMbCVG/XmeF8ElIjX5PnB7586dqn/v27cPtra2qts5OTk4cOAAvLy8tFpcYVq8eDFmzpyJ+/fvo2rVqpg3bx6aNGmi67KIqAhJJBLM+LgGguYewfk7T7DsyE182dhT12URkZ7Id0jq1KkTgBc/VHr16qV2n0wmg5eXF2bPnq3V4grL5s2bMWLECCxevBiNGjXCsmXL0KZNG0RFRaFMmTK6Lo+IipCrrQUmdqiK0VvOY97+q2ha3k7XJRGRnsj37jalUgmlUokyZcogMTFRdVupVCIzMxMxMTFo3759YdaqNXPmzEG/fv3Qv39/VKlSBfPmzYOHhweWLFmi69KISAc61y6NgCrOUOQIjP3rErjXjYiAdzhPUmxsbGHUUWSysrIQGRmJcePGqY0HBQUhLCxM42MyMzPVrk+XmpoK4MUJNRUKhdZqy12WNpepT4y9P8D4ezTm/iZ3qIzIW48RnZCGYFMpWhthj4Bxv4aA8fcHGH+PhdlfQZcpEfk4L/+vv/6KL7/8Eubm5vj111/fOHfYsGEFKqCoxcfHo3Tp0jh+/DgaNmyoGp8yZQrWrFmDmJiYPI+ZOHEiJk2alGd8w4YNsLTkSeiIjMWZRxKsuWYCqURgVLUceJTQdUVEpE3p6eno0aMHUlJSYGNj89b5+QpJ3t7eOH36NOzt7eHt7f36hUkkuHnzZsEqLmK5ISksLAwNGjRQjf/yyy9Yu3Ytrly5kucxmrYkeXh44NGjR/l6kvNLoVAgJCQEgYGBrz0HlSEz9v4A4+/R2PsTQmDIxnMIjn6I8o5W2DGoAeSm73z1Jr1k7K+hsfcHGH+PhdlfamoqHBwc8h2S8rW77eVdbIa+u83BwQEmJiZISEhQG09MTISzs7PGx8jlcsjl8jzjMpmsUN6ghbVcfWHs/QHG36Mx9ze5Y1WEXT+E6w+fYdHhWHzbuvLbH2SAjPk1BIy/P8D4eyyM/gq6POP6L1I+mJmZoU6dOggJCVEbDwkJUdv9RkTFk72VGbp6vzhye9nhGzhzO1nHFRGRrhQ4JH388ceYNm1anvGZM2fik08+0UpRhW3UqFH4/fffsXLlSkRHR2PkyJG4ffs2vvrqK12XRkR6oKa9QMearlAKYPSf55GRlfcEukRk/Aockg4fPox27drlGW/dujWOHDmilaIKW7du3TBv3jxMnjwZtWrVwpEjR7B79254evIkckT0wo/tKsPFxhyxj55h+t68xyoSkfErcEh6+vQpzMzM8ozLZDLVV+MNwaBBgxAXF4fMzExERkaiadOmui6JiPSIrYUM0z+uAQBYHRaH49cf6bgiIipqBQ5J1apVw+bNm/OMb9q0CT4+PlopiohIH/hXdETPD16chX/MlvNIyTDO89IQkWYFPpnkjz/+iC5duuDGjRto0aIFAODAgQPYuHEjtmzZovUCiYh06bu2VXDs+iPcSkrHpJ2XMadbLV2XRERFpMBbkj788EPs2LED169fx6BBgzB69GjcvXsX+/fvV13fjYjIWFjJTTGna01IJcC2s/ew5+J9XZdEREWkwFuSAKBdu3YaD94mIjJGdTzt8JV/OSwOvYHvtl9EHc9ScLIx13VZRFTI3vk8SZGRkVi3bh3Wr1+Ps2fParMmIiK9MyKgInxcbZCcrsDYvy4gHxcrICIDV+CQlJiYiBYtWqBu3boYNmwYhgwZgjp16qBly5Z4+PBhYdRIRKRzZqZSzOteC2amUoTGPMS6k7d1XRIRFbICh6ShQ4ciNTUVly9fxuPHj5GcnIxLly4hNTVV7y9uS0T0Pio6W2Pcf5cp+WVXFK4nPtVxRURUmAockvbu3YslS5agSpUqqjEfHx8sWrQIe/bs0WpxRET6pndDLzQu74DnCiVGbj6HrGylrksiokJS4JCkVCo1XiBOJpNBqeQPCyIyblKpBLM+qQlbCxku3kvB/ANXdV0SERWSAoekFi1aYPjw4YiPj1eN3bt3DyNHjkTLli21WhwRkT5ysTXHlI+qAwAWh97AyZtJOq6IiApDgUPSwoULkZaWBi8vL5QrVw7ly5eHt7c30tLSsGDBgsKokYhI77Sr4YqP67hDCGDk5nNISefZuImMTYHPk+Th4YEzZ84gJCQEV65cgRACPj4+CAgIKIz6iIj01sQPqyIi7jFuJaXjux0XsfBTX0gkEl2XRURa8k4nkwSAwMBABAYGarMWIiKDUkJuivndffHxkjDsunAf/hUd0dXPQ9dlEZGW5Csk/frrr/leIE8DQETFSS2PkhgZWBEz98Vg4s7LqOtlB28HK12XRURakK+QNHfu3HwtTCKRMCQRUbHzlX85HLn6ECdjH2PEprPY+nVDyEze+YIGRKQn8hWSYmNjC7sOIiKDZSKVYG63Wmgz/yjO303B3JCrGPvfSSeJyHC98391srKyEBMTg+zsbG3WQ0RkkNxKWmBq5xenBVhy+AbCbjzScUVE9L4KHJLS09PRr18/WFpaomrVqrh9+8X1i4YNG4Zp06ZpvUAiIkPRtroruvl5qE4LkPQ0U9clEdF7KHBIGj9+PM6fP4/Q0FCYm5urxgMCArB582atFkdEZGgmfOiDco5WeJCaiTFbzkMIoeuSiOgdFTgk7dixAwsXLkTjxo3Vzgfi4+ODGzduaLU4IiJDY2lmioU9asPMVIpDMQ+x4hiP6SQyVAUOSQ8fPoSTk1Oe8WfPnvEkakREAKq42uDH9j4AgOl7r+DC3Se6LYiI3kmBQ1LdunWxa9cu1e3cYLR8+XI0aNBAe5URERmwzz4og9ZVXaDIERi68SzSnvOyJUSGpsBn3J46dSpat26NqKgoZGdnY/78+bh8+TLCw8Nx+PDhwqiRiMjgSCQSTO9SAxfvpeBWUjq+334J87vX4hZ3IgOS7y1J586dAwA0bNgQx48fR3p6OsqVK4fg4GA4OzsjPDwcderUKaw6iYgMjq2lDL9+WgsmUgl2no/HltN3dV0SERVAvrck1a5dG76+vujfvz969OiBNWvWFGZdRERGoY6nHUYHVcSMvTH4385LqOlREpVcrHVdFhHlQ763JB0/fhy1a9fGuHHj4Orqis8//xyHDh0qzNqIiIzCV03LoUkFBzxXKPH1+kg8zeRJeIkMQb5DUoMGDbB8+XIkJCRgyZIluHPnDgICAlCuXDn88ssvuHuXm5GJiDSRSiWY160WXG3NcfPhM3z71wWeP4nIABT4220WFhbo1asXQkNDcfXqVXz66adYtmwZvL290bZt28KokYjI4NmXkGNhj9owlUqw68J9rAmL03VJRPQW73WZ6nLlymHcuHH4/vvvYWNjg3379mmrLiIio1PHsxS+a1sFAPDL7micuZ2s44qI6E3eOSQdPnwYvXr1gouLC8aOHYvOnTvj+PHj2qyNiMjo9GnkhXbVXaHIERi8/gweP8vSdUlE9BoFCkl37tzBTz/9hHLlyqF58+a4ceMGFixYgPj4eCxfvhz169cvrDqJiIyCRCLBtC7VUdbBCvdTnmP4prPIUfL4JCJ9lO+QFBgYCG9vbyxevBgff/wxoqOjcezYMfTp0wdWVlaFWSMRkVGxNpdh8We1YS6T4ui1R1hw8JquSyIiDfIdkiwsLPDXX3/h7t27mD59OipVqlSYdRERGbXKLjb4pVN1AMD8A9dw5OpDHVdERK/Kd0jauXMnOnbsCBMTk8Ksh4io2OhSxx2f1isDIYDhm84i/kmGrksiope817fbiIjo/Uzo4IOqbjZITldg8IYzyMpW6rokIvoPQxIRkQ6Zy0ywpGcd2Jib4uztJ5i6J1rXJRHRfxiSiIh0rIy9JWZ3rQUAWHU8Drsu3NdtQUQEgCGJiEgvBPo44yv/cgCAsVvP48bDpzquiIgYkoiI9MSYoIr4wNsOz7Jy8PW6SKRn8UK4RLrEkEREpCdMTaRY0MMXjtZyXH3wFGO38kK4RLrEkEREpEecrM2x6L8L4f574T4Wh97QdUlExRZDEhGRnqnnbYdJHasCAGYFx+BA9AMdV0RUPDEkERHpoZ4feKLnB7knmjyH64lpui6JqNhhSCIi0lMTOlRFPW87PM3MxoA/IpGSodB1SUTFCkMSEZGeMjOVYnHP2ihd0gKxj55h2MazyFHyQG6iosKQRESkxxxKyLHs8zowl0lx+OpDzNh7RdclERUbDElERHquWmlbzPy4JgBg2ZGb2HH2no4rIioeGJKIiAxAh5puGNTsxRm5v/3rAi7cfaLbgoiKAYYkIiIDMSaoElpWdkJmthJf/hGJxLTnui6JyKgxJBERGQipVIK53WuhnKMVElKf4+t1Z5CZnaPrsoiMFkMSEZEBsTGXYfkXfrA2N0XkrWT8b8dlXrqEqJAwJBERGZiyjiWw4FNfSCXA5tN3sPbELV2XRGSUGJKIiAxQs0pO+LZ1ZQDApH+iEHbjkY4rIjI+DElERAbqy6Zl0amWG3KUAl+tjeSlS4i0jCGJiMhASSQSTOtSA7XLlETq82z0WR2BR08zdV0WkdFgSCIiMmDmMhMs/8IPZewscedxBvqvOY2MLH7jjUgbGJKIiAycfQk5Vvepi5KWMpy78wQjNvMab0TawJBERGQEyjqWwG+f+8HMRIp9lx9g6u5oXZdEZPAYkoiIjEQ9bzvM/KQGAOD3Y7H4IzxOtwURGTiGJCIiI9KxVmmMCaoIAJi48zIORD/QcUVEhoshiYjIyAxuXh5d/dyhFMDQjWdx6V6KrksiMkgMSURERkYikeCXj6qjcXkHpGfloO/qCNx7kqHrsogMDkMSEZERkplIsfiz2qjoXAKJaZnouyoCqc8Vui6LyKAwJBERGSkbcxlW9akHR2s5Yh6kYfD6M1DkKHVdFpHBYEgiIjJipUtaYGWvurCQmeDotUeY8E80BE+hRJQvDElEREauurstFnzqC6kE2BJ5DyH3JLouicggGERIiouLQ79+/eDt7Q0LCwuUK1cOEyZMQFZWltq827dvo0OHDrCysoKDgwOGDRuWZ87Fixfh7+8PCwsLlC5dGpMnT4bgf6uIyMgF+DhjQoeqAIBdd0zw5+m7Oq6ISP+Z6rqA/Lhy5QqUSiWWLVuG8uXL49KlSxgwYACePXuGWbNmAQBycnLQrl07ODo64tixY0hKSkKvXr0ghMCCBQsAAKmpqQgMDETz5s0RERGBq1evonfv3rCyssLo0aN12SIRUaHr1dAL95Kf4bejcfhxZxTsSpijTXVXXZdFpLcMIiS1bt0arVu3Vt0uW7YsYmJisGTJElVICg4ORlRUFO7cuQM3NzcAwOzZs9G7d2/88ssvsLGxwfr16/H8+XOsXr0acrkc1apVw9WrVzFnzhyMGjUKEgk3QRORcRsTWAEXY24iPFGK4ZvOwdpchsYVHHRdFpFeMoiQpElKSgrs7OxUt8PDw1GtWjVVQAKAVq1aITMzE5GRkWjevDnCw8Ph7+8PuVyuNmf8+PGIi4uDt7e3xnVlZmYiMzNTdTs1NRUAoFAooFBo7yu1ucvS5jL1ibH3Bxh/j8beH2D8PWZnZ6NrWSWs7Z0RHP0QX649jTW966CWR0ldl6YVxv76AcbfY2H2V9BlGmRIunHjBhYsWIDZs2erxhISEuDs7Kw2r1SpUjAzM0NCQoJqjpeXl9qc3MckJCS8NiRNnToVkyZNyjMeHBwMS0vL92lFo5CQEK0vU58Ye3+A8fdo7P0Bxt2jVAK0srmPW7ZSxKQAvVaexLCqOXDV/o8znTHm1y+XsfdYGP2lp6cXaL5OQ9LEiRM1ho+XRUREwM/PT3U7Pj4erVu3xieffIL+/furzdW0u0wIoTb+6pzcg7bftKtt/PjxGDVqlOp2amoqPDw8EBQUBBsbmzfWXxAKhQIhISEIDAyETCbT2nL1hbH3Bxh/j8beH2D8Peb216ZVIFoGStBrdSTO303BqptW2DSgHtxLWei6xPdi7K8fYPw9FmZ/uXuC8kunIWnIkCHo3r37G+e8vOUnPj4ezZs3R4MGDfDbb7+pzXNxccHJkyfVxpKTk6FQKFRbi1xcXFRblXIlJiYCQJ6tUC+Ty+Vqu+hyyWSyQnmDFtZy9YWx9wcYf4/G3h9g/D3KZDJYymRY3aceui4Lx7XEp+izJhJbvmoIR+u8P+8MjbG/foDx91gY/RV0eTo9BYCDgwMqV678xj/m5uYAgHv37qFZs2aoXbs2Vq1aBalUvfQGDRrg0qVLuH//vmosODgYcrkcderUUc05cuSI2mkBgoOD4ebmlmc3HBFRcVDKygxr+30A91IWiEtKR6+Vp3j5EqL/GMR5kuLj49GsWTN4eHhg1qxZePjwIRISEtS2CgUFBcHHxweff/45zp49iwMHDmDMmDEYMGCAapdYjx49IJfL0bt3b1y6dAnbt2/HlClT+M02IirWXGzNsbbfB3AoYYao+6nov/o0MrJydF0Wkc4ZREgKDg7G9evXcfDgQbi7u8PV1VX1J5eJiQl27doFc3NzNGrUCF27dkWnTp1UpwgAAFtbW4SEhODu3bvw8/PDoEGDMGrUKLXjjYiIiiNvByus6VsP1nJTnIp7jMEbeJ03IoP4dlvv3r3Ru3fvt84rU6YM/v333zfOqV69Oo4cOaKlyoiIjEdVN1us6F0Xn684iYNXEvHNlvOY07UWpFJuaafiySC2JBERUdGo522HJZ/VhqlUgh3n4jH53yheuomKLYYkIiJS06KyM2Z9UhMAsDosDtP2XGFQomKJIYmIiPLo5FsaP3WqBgBYduQmpu1lUKLihyGJiIg0+ry+JyZ3rAoAWHb4Jmbsi2FQomKFIYmIiF7riwZemNjBBwCwJPQGZgUzKFHxwZBERERv1LuRNyb8F5QWHbqB2cFXGZSoWGBIIiKit+rTyBs/tn8RlBYeuo65+6/puCKiwseQRERE+dKvsTd+aFcFAPDrgWuYG3JVxxURFS6GJCIiyrf+TcqqgtL8A9cwbz+DEhkvhiQiIiqQ/k3K4ru2lQEA8/Zfw68HuOuNjBNDEhERFdiXTcthXJsXQWlOyFUsPMigRMaHIYmIiN7JV/7l8G3rF0FpVvBVLDp0XccVEWkXQxIREb2zr5uVwzetKgEAZu6Lwbz9PD0AGQ+GJCIiei+Dm5dXBaV5+69h8r9RUCoZlMjwMSQREdF7G9y8vOqEk6uOx+GbrReQnaPUcVVE74chiYiItKJPI2/M6VoTJlIJ/jpzF1+vP4Pnihxdl0X0zhiSiIhIazrXdsfSz+rAzFSKkKgH6Ls6Ak8zs3VdFtE7YUgiIiKtCvRxxuo+dWFlZoKwG0noufwEkp9l6bosogJjSCIiIq1rWM4BG7+sj1KWMpy/m4Kuy8KRkPJc12URFQhDEhERFYoa7iWx5asGcLExx7XEp/h4aRjiHj3TdVlE+caQREREhaa8kzW2ft0AXvaWuJucgY+XhiP6fqquyyLKF4YkIiIqVO6lLLHlq4ao4mqDR08z0W1ZOCJvPdZ1WURvxZBERESFztFajk1f1oefZymkPs/GZ7+fQmhMoq7LInojhiQiIioSthYy/NGvHvwrOiJDkYN+a05jTVicrssiei2GJCIiKjKWZqZY/oUfutR2R45SYMLOy/hxxyUoeHZu0kMMSUREVKTMTKWY9UkNjGtTGRIJsPbELfRZFYGUdIWuSyNSw5BERERFTiKR4Cv/clj6WR1Ympng2PVH+GjJccTyFAGkRxiSiIhIZ1pVdcGWrxrAzdYcNx8+Q6dFxxF+I0nXZREBYEgiIiIdq+pmix1DGqGmR0mkZCjw+YqT2HTqtq7LImJIIiIi3XOyNsfmL+ujQ003ZCsFxm27iJ/+jUKOUui6NCrGGJKIiEgvmMtM8Gv3WhgZUBEAsOJYLAb8cRppz3lAN+kGQxIREekNiUSC4QEVsLCHL+SmUhy8koiPl4TjzuN0XZdGxRBDEhER6Z32NdyweWADOFrLEfMgDZ0WHcfpOF7KhIoWQxIREemlWh4lsXNII1R1s0HSsyz0WH4S287c1XVZVIwwJBERkd5ytbXAlq8aoFVVZ2TlKDHqz/OYsfcKlDygm4oAQxIREek1SzNTLOlZB4OalQMALA69ga/XRyI9K1vHlZGxY0giIiK9J5VKMLZ1ZczpWhNmJlLsu/wAnywNR/yTDF2XRkaMIYmIiAxG59ru2DDgA9hbmeFyfCra/noU+y4n6LosMlIMSUREZFD8vOywY3AjVC9tiyfpCgxcG4kfd1zCc0WOrksjI8OQREREBsfDzhJ/fd0QXzYtCwBYe+IWOi48jmsPnuq4MjImDElERGSQzEyl+K5tFazpWw8OJcwQ8yANHy09geMPJBCC336j98eQREREBs2/oiP2DG+KphUdkZmtxJ83TTBk03k8Sc/SdWlk4BiSiIjI4Dlay7G6d12Ma10RJhKB4KhEtJ1/FKdieZZuencMSUREZBSkUgn6NfLCiGo58LSzRHzKc3T/LRzz9l9Fdo5S1+WRAWJIIiIio1KmBLBjUH10rl0aSgHM238NPZaf5DmVqMAYkoiIyOiUkJtiTtdamNutJqzMTHAq7jHazD+KvZd4TiXKP4YkIiIyWh/5umP38Cao6W6LlAwFvloXie+3X+Q5lShfGJKIiMioedpbYctXDTHQ/8U5ldafvI0PFx5DTEKajisjfceQRERERs/MVIrxbapgbb96cCghx9UHT/HhwmNYe+IWz6lEr8WQRERExUaTCo7YO6IJmlV6cU6lH3dcwhcrTyH20TNdl0Z6iCGJiIiKFYcScqzsVRc/tKsCM1Mpjl57hFbzjmDe/qs8VonUMCQREVGxI5VK0L9JWQSPaIomFRyQla3EvP3X0HreERy99lDX5ZGeYEgiIqJiy8vBCn/0rYeFPXzhZC1HXFI6Pl9xCkM3nkVi6nNdl0c6xpBERETFmkQiQfsabjgw2h99GnlBKgH+OR+PlrMPY01YHHKUPLC7uGJIIiIiAmBtLsOEDlWxc0hj1HS3RVpmNibsvIxOi47jwt0nui6PdIAhiYiI6CXVStti26BG+KlTNVibm+LivRR0XHQc//v7ElKfK3RdHhUhhiQiIqJXmEgl+Ly+Jw6OboaPfEtDCOCP8FtoOfsw/j53j+dWKiYYkoiIiF7D0VqOud1qYUP/D1DWwQoP0zIxfNM5fLbiJG4+fKrr8qiQMSQRERG9RcPyDtgzoglGB1aE3FSK49eT0HreUcwJ4bmVjBlDEhERUT7ITU0wtGUFBI9sCv+KjsjKUeLXA9fQat4R/HshHkp+C87oMCQREREVgKe9FVb3qYvFPWvD2UaOW0npGLLhLFrPP4JdF+4zLBkRhiQiIqICkkgkaFvdFftH+WNEQAVYm5vi6oOnGLzhDNr+ehR7LjIsGQOGJCIiondkbS7DiICKOPZtCwxvWQHWclNcSUjD1+tfhKW9lxiWDBlDEhER0XuytZBhZOCLsDSsRXmU+C8sfbXuDNotOIZ9lxN42gADxJBERESkJbaWMowKqoRj3zbH0P/CUvT9VAxcG4l2vx5DMMOSQWFIIiIi0rKSlmYYHVQJR8c2x+Dm5WBlZoKo+6n4cm0kOiw8hv1RDxiWDABDEhERUSEpZWWGb1pVxrFvW2BQsxdh6dK9VPT/4zQ+XHgcB6IZlvSZwYWkzMxM1KpVCxKJBOfOnVO77/bt2+jQoQOsrKzg4OCAYcOGISsrS23OxYsX4e/vDwsLC5QuXRqTJ0/mG5SIiApVKSszjG1dGUe/bYGvm5WDpZkJLt5LQb81p9Fp0XEcvPKAB3jrIVNdF1BQY8eOhZubG86fP682npOTg3bt2sHR0RHHjh1DUlISevXqBSEEFixYAABITU1FYGAgmjdvjoiICFy9ehW9e/eGlZUVRo8erYt2iIioGLGzMsO3rSujf2Nv/Hb0Jv4Iu4Xzd1PQd/VplLGzRLe6Hvi4jjucbcx1XSrBwLYk7dmzB8HBwZg1a1ae+4KDgxEVFYV169bB19cXAQEBmD17NpYvX47U1FQAwPr16/H8+XOsXr0a1apVQ+fOnfHdd99hzpw53JpERERFxr6EHOPbVMGxb5tjYNOysJab4vbjdMzcF4MGUw+g/5oIhEQ9QHaOUtelFmsGsyXpwYMHGDBgAHbs2AFLS8s894eHh6NatWpwc3NTjbVq1QqZmZmIjIxE8+bNER4eDn9/f8jlcrU548ePR1xcHLy9vTWuOzMzE5mZmarbuaFLoVBAoVBoq0XVsrS5TH1i7P0Bxt+jsfcHGH+P7E+/2MilGBNYHoP8vbD38gNsibyH07eeYH90IvZHJ8LJWo4uvm7oUqc0PO1e/O4ztB4LqjD7K+gyJcIANqEIIdC2bVs0atQIP/zwgyrQnD17FrVq1QIAfPnll4iLi0NwcLDaY+VyOVavXo1PP/0UQUFB8PLywm+//aa6Pz4+HqVLl0ZYWBgaNGigcf0TJ07EpEmT8oxv2LBBY2AjIiJ6Vw8ygBMPpDj1UIKn2RLVeAUbJRo4C9SwE5AZ1H4g/ZGeno4ePXogJSUFNjY2b52v0y1JrwsfL4uIiEBYWBhSU1Mxfvz4N86VSCR5xoQQauOvzsnNiJoem2v8+PEYNWqU6nZqaio8PDwQFBSUryc5vxQKBUJCQhAYGAiZTKa15eoLY+8PMP4ejb0/wPh7ZH+GoQ+ArGwlDsY8xJ+n7+LYjSRcS5XiWipgY26KWiWzMKrjB6jqXkrXpWpdYb6GuXuC8kunIWnIkCHo3r37G+d4eXnh559/xokTJ9R2kwGAn58fevbsiTVr1sDFxQUnT55Uuz85ORkKhQLOzs4AABcXFyQkJKjNSUxMBADVHE3kcnmedQOATCYrlA9hYS1XXxh7f4Dx92js/QHG3yP7038yGdChljs61HLH3eR0bDl9F1tO30F8ynMcSZDiyLII1PIoie51PdC+phtKyA3mCJp8KYzXsKDL0+kz6uDgAAcHh7fO+/XXX/Hzzz+rbsfHx6NVq1bYvHkzPvjgAwBAgwYN8Msvv+D+/ftwdXUF8OJgbrlcjjp16qjmfPfdd8jKyoKZmZlqjpubG7y8vLTcHRERkXa4l7LEyMCKGNayAkKvJGD+v6cR9cQE5+48wbk7TzD53yi0re6KgCpOaFjeATbmhh0Q9YVBxM4yZcqo3S5RogQAoFy5cnB3dwcABAUFwcfHB59//jlmzpyJx48fY8yYMRgwYIBql1iPHj0wadIk9O7dG9999x2uXbuGKVOm4H//+98bd7cRERHpAxOpBE0rOOBpJSXqNW2OnRceYHPEHdx89AxbI+9ia+RdmEglqF2mJPwrOqJpRUdUc7OFVMrfce/CIEJSfpiYmGDXrl0YNGgQGjVqBAsLC/To0UPtdAG2trYICQnB4MGD4efnh1KlSmHUqFFqxxsREREZAocScgz0L4cvm5bFqdjH2HMpAUeuPsTNR88QEZeMiLhkzAq+CjsrMzSp4AD/io5oUsERjtZ5Dx8hzQwyJHl5eWk8r1GZMmXw77//vvGx1atXx5EjRwqrNCIioiIlkUjwQVl7fFDWHgBw53E6Dl99iCNXHyLsRhIeP8vC3+fi8fe5eABAVTcb1VamOp6lIDPhV+VexyBDEhEREWnmYWeJz+p74rP6nlDkKHHmVvKL0HTtIS7dS8Xl+Bd/FofeQAm5KRqUs4d/RUf4V3SEhx1Pa/MyhiQiIiIjJTORqrYyjW1dGQ/TMnHs+kMcjnmIo9ceIelZFkKiHiAk6gEAwMPOApWcbVDRuQQqOlujgnMJlHMsAXOZiY470Q2GJCIiomLC0VqOj3zd8ZGvO5RKgcvxqThy7UVoirydjDuPM3DncQb2Rz9QPUYqAcrYWaKCs/X/hycna5R1tDL68MSQREREVAxJpRJUd7dFdXdbDG5eHqnPFbh0LwXXHjzF1QdpL/5OTMOTdAXiktIRl5Su2uIEvAhPXvZWqKDa6vQiRHmUsoSlmYlRfGucIYmIiIhgYy5Dw3IOaFju/89fKITAo6dZuPYgDVcfpOFq4tP//v0UKRkK3Hz0DDcfPcO+yw/UliWVACXkprA2l8Ha3PS/f5uihLkMJeSmsPlvrITqvv+fZ2EKpGYBSqXur5rGkEREREQaSSQSOFrL4WgtR8Py6uHpYVomruZudUpMU/077Xk2lAJIfZ6N1OfZ77hmU7RvkwMNF7soUgxJREREVCASiQRONuZwsjFH4wrq4SlDkYOnz7ORlpmNtOfZePo8G08zFUhV/fvFn7Tnihf3q83LRupzBZ5lKmBlpvvjnRiSiIiISCskEgkszUxhaWYKp3dchkKhwK5du/XimCaeQYqIiIj0ih7kIwAMSUREREQaMSQRERERacCQRERERKQBQxIRERGRBgxJRERERBowJBERERFpwJBEREREpAFDEhEREZEGDElEREREGjAkEREREWnAkERERESkAUMSERERkQYMSUREREQamOq6AEMkhAAApKamanW5CoUC6enpSE1NhUwm0+qy9YGx9wcYf4/G3h9g/D2yP8Nn7D0WZn+5v7dzf4+/DUPSO0hLSwMAeHh46LgSIiIiKqi0tDTY2tq+dZ5E5DdOkYpSqUR8fDysra0hkUi0ttzU1FR4eHjgzp07sLGx0dpy9YWx9wcYf4/G3h9g/D2yP8Nn7D0WZn9CCKSlpcHNzQ1S6duPOOKWpHcglUrh7u5eaMu3sbExyjd+LmPvDzD+Ho29P8D4e2R/hs/Yeyys/vKzBSkXD9wmIiIi0oAhiYiIiEgDhiQ9IpfLMWHCBMjlcl2XUiiMvT/A+Hs09v4A4++R/Rk+Y+9Rn/rjgdtEREREGnBLEhEREZEGDElEREREGjAkEREREWnAkERERESkAUPSe1i8eDG8vb1hbm6OOnXq4OjRo2+cf/jwYdSpUwfm5uYoW7Ysli5dmmfOX3/9BR8fH8jlcvj4+GD79u0FXq8QAhMnToSbmxssLCzQrFkzXL582SD6mzp1KurWrQtra2s4OTmhU6dOiImJUZvTu3dvSCQStT/169c3iP4mTpyYp3YXFxe1Odp6/XTVo5eXV54eJRIJBg8erJqjr6/h5cuX0aVLF1UP8+bNe6f16utnMD/9FeVnUFc9FuXnUBf9GfJncPny5WjSpAlKlSqFUqVKISAgAKdOnSrwerX2c1TQO9m0aZOQyWRi+fLlIioqSgwfPlxYWVmJW7duaZx/8+ZNYWlpKYYPHy6ioqLE8uXLhUwmE1u3blXNCQsLEyYmJmLKlCkiOjpaTJkyRZiamooTJ04UaL3Tpk0T1tbW4q+//hIXL14U3bp1E66uriI1NVXv+2vVqpVYtWqVuHTpkjh37pxo166dKFOmjHj69KlqTq9evUTr1q3F/fv3VX+SkpLy3Zsu+5swYYKoWrWqWu2JiYlq69LG66fLHhMTE9X6CwkJEQDEoUOHVHP09TU8deqUGDNmjNi4caNwcXERc+fOfaf16utnMD/9FdVnUJc9FtXnUFf9GfJnsEePHmLRokXi7NmzIjo6WvTp00fY2tqKu3fvFmi92vo5ypD0jurVqye++uortbHKlSuLcePGaZw/duxYUblyZbWxgQMHivr166tud+3aVbRu3VptTqtWrUT37t3zvV6lUilcXFzEtGnTVPc/f/5c2NraiqVLl+p9f69KTEwUAMThw4dVY7169RIdO3bMbysa6aq/CRMmiJo1a762Lm29fkLoz2s4fPhwUa5cOaFUKlVj+voavszT01PjLyBD/gy+7HX9vaqwPoNC6K7Hovoc6straKifQSGEyM7OFtbW1mLNmjX5Xq82f45yd9s7yMrKQmRkJIKCgtTGg4KCEBYWpvEx4eHheea3atUKp0+fhkKheOOc3GXmZ72xsbFISEhQmyOXy+Hv7//a2vSlP01SUlIAAHZ2dmrjoaGhcHJyQsWKFTFgwAAkJibmqzdA9/1du3YNbm5u8Pb2Rvfu3XHz5k3Vfdp4/fShx5frWLduHfr27ZvnYtD6+BpqY736/Bl8F4XxGQR032Nhfw513d/LdRjyZzA9PR0KhUL1/iuqz2AuhqR38OjRI+Tk5MDZ2Vlt3NnZGQkJCRofk5CQoHF+dnY2Hj169MY5ucvMz3pz/y5IbfrS36uEEBg1ahQaN26MatWqqcbbtGmD9evX4+DBg5g9ezYiIiLQokULZGZm6n1/H3zwAf744w/s27cPy5cvR0JCAho2bIikpCTVMnIfl9/a9K3Hl+3YsQNPnjxB79691cb19TXUxnr1+TNYUIX1GQR022NRfA715TU09M/guHHjULp0aQQEBOR7vdr6OQoApgWaTWpeTeVCiDxjb5v/6nh+lqmtOW+jq/5yDRkyBBcuXMCxY8fUxrt166b6d7Vq1eDn5wdPT0/s2rULnTt3fkNHb6+3sPtr06aN6t/Vq1dHgwYNUK5cOaxZswajRo1659oKUnNRvoYrVqxAmzZt4Obmpjauz6+httarr5/Bgijsz+Drai7sHovyc6jr19CQP4MzZszAxo0bERoaCnNz8wKvVxuvH7ckvQMHBweYmJjkSaSJiYl5kmsuFxcXjfNNTU1hb2//xjm5y8zPenO/oVGQ2vSlv5cNHToUO3fuxKFDh+Du7v7Gel1dXeHp6Ylr1669tTdAP/rLZWVlherVq6tq18brB+hHj7du3cL+/fvRv3//t9arL6+hNtarz5/BgijMzyCgHz3mKozPoT70Z8ifwVmzZmHKlCkIDg5GjRo1CrRebf0cBRiS3omZmRnq1KmDkJAQtfGQkBA0bNhQ42MaNGiQZ35wcDD8/Pwgk8neOCd3mflZr7e3N1xcXNTmZGVl4fDhw6+tTV/6A14k/SFDhmDbtm04ePAgvL2931pvUlIS7ty5A1dXV73v71WZmZmIjo5W1a6N109fely1ahWcnJzQrl27t9arL6+hNtarz5/B/CiKzyCg2x5fVRifQ33oz1A/gzNnzsRPP/2EvXv3ws/Pr8Dr1dbPUQA8BcC7yv0K4ooVK0RUVJQYMWKEsLKyEnFxcUIIIcaNGyc+//xz1fzcrz6OHDlSREVFiRUrVuT56uPx48eFiYmJmDZtmoiOjhbTpk177SkAXrdeIV589dHW1lZs27ZNXLx4UXz66afv/NXVou7v66+/Fra2tiI0NFTtq6np6elCCCHS0tLE6NGjRVhYmIiNjRWHDh0SDRo0EKVLlzaI/kaPHi1CQ0PFzZs3xYkTJ0T79u2FtbW11l8/XfYohBA5OTmiTJky4ttvv81Tlz6/hpmZmeLs2bPi7NmzwtXVVYwZM0acPXtWXLt2Ld/rFUJ/P4P56a+oPoO67LGoPoe66k8Iw/0MTp8+XZiZmYmtW7eqvf/S0tLyvV4htPdzlCHpPSxatEh4enoKMzMzUbt27TxfkfX391ebHxoaKnx9fYWZmZnw8vISS5YsybPMLVu2iEqVKgmZTCYqV64s/vrrrwKtV4gXX3+cMGGCcHFxEXK5XDRt2lRcvHjRIPoDoPHPqlWrhBBCpKeni6CgIOHo6ChkMpkoU6aM6NWrl7h9+7ZB9Jd7rg6ZTCbc3NxE586dxeXLl9XmaOv101WPQgixb98+AUDExMTkuU+fX8PY2FiN779Xl2Oon8H89FeUn0Fd9ViUn0NdvUcN9TPo6empsb8JEybke71CaO/1kwjx31FTRERERKTCY5KIiIiINGBIIiIiItKAIYmIiIhIA4YkIiIiIg0YkoiIiIg0YEgiIiIi0oAhiYiIiEgDhiQiKnKhoaGQSCR48uTJa+esXr0aJUuWLLKaNOnduzc6der0xjn56UVbkpKS4OTkhLi4uHzNz8zMRJkyZRAZGVm4hREZKYYkInonCQkJGD58OMqXLw9zc3M4OzujcePGWLp0KdLT09/42IYNG+L+/fuwtbUtomrfzfz587F69WrV7WbNmmHEiBFqc4qyl6lTp6JDhw7w8vLK13y5XI4xY8bg22+/LdzCiIyUqa4LICLDc/PmTTRq1AglS5bElClTUL16dWRnZ+Pq1atYuXIl3Nzc8OGHH2p8rEKhgJmZmepK3fosP8GnqHrJyMjAihUrsHv37gI9rmfPnvjmm28QHR2NKlWqFFJ1RMaJW5KIqMAGDRoEU1NTnD59Gl27dkWVKlVQvXp1dOnSBbt27UKHDh1UcyUSCZYuXYqOHTvCysoKP//8s8ZdVKtXr0aZMmVgaWmJjz76CElJSW+sIS4uDhKJBJs2bULDhg1hbm6OqlWrIjQ0VG3e4cOHUa9ePcjlcri6umLcuHHIzs5W3b9161ZUr14dFhYWsLe3R0BAAJ49ewZAfXdb7969cfjwYcyfPx8SiQQSiQRxcXEae/nrr79QtWpVyOVyeHl5Yfbs2Wo1eXl5YcqUKejbty+sra1RpkwZ/Pbbb2/sd8+ePTA1NUWDBg1UY5MnT4abm5vac/Xhhx+iadOmUCqVAAB7e3s0bNgQGzdufOPyiUiDAl/tjYiKtUePHgmJRCKmTp2ar/kAhJOTk1ixYoW4ceOGiIuLE4cOHRIARHJyshBCiBMnTqiWGRMTI+bPny9KliwpbG1tX7vc3At9uru7i61bt4qoqCjRv39/YW1tLR49eiSEEOLu3bvC0tJSDBo0SERHR4vt27cLBwcH1cUy4+PjhampqZgzZ46IjY0VFy5cEIsWLVJdcbxXr16iY8eOQgghnjx5Iho0aCAGDBigujJ5dnZ2nl5Onz4tpFKpmDx5soiJiRGrVq0SFhYWqgvECvHiIp52dnZi0aJF4tq1a2Lq1KlCKpWK6Ojo1/Y7fPhw0bp1a7Wx7Oxs0aBBA9GpUychhBBLliwRtra2aldDF0KIsWPHimbNmr3pZSIiDRiSiKhATpw4IQCIbdu2qY3b29sLKysrYWVlJcaOHasaByBGjBihNvfVYPHpp5/mCQDdunXLV0iaNm2aakyhUAh3d3cxffp0IYQQ3333nahUqZJQKpWqOYsWLRIlSpQQOTk5IjIyUgDIEypyvRyShBDC399fDB8+/I299OjRQwQGBqrN+eabb4SPj4/qtqenp/jss89Ut5VKpXBycspzRfSXdezYUfTt2zfP+I0bN4S1tbX49ttvhaWlpVi3bl2eOfPnzxdeXl6vXTYRacbdbUT0TiQSidrtU6dO4dy5c6hatSoyMzPV7vPz83vjsqKjo9V2IwHIc/t1Xp5namoKPz8/REdHqy335VobNWqEp0+f4u7du6hZsyZatmyJ6tWr45NPPsHy5cuRnJycr/W+qZdGjRqpjTVq1AjXrl1DTk6OaqxGjRqqf0skEri4uCAxMfG1y83IyIC5uXme8bJly2LWrFmYPn06OnTogJ49e+aZY2Fh8daD6YkoL4YkIiqQ8uXLQyKR4MqVK2rjZcuWRfny5WFhYZHnMVZWVm9cphBCqzXmhiIhRJ4wl7suiUQCExMThISEYM+ePfDx8cGCBQtQqVIlxMbGvvO637TOl8lksjw15x5HpImDg8NrA9yRI0dgYmKCuLg4teOtcj1+/BiOjo75KZ+IXsKQREQFYm9vj8DAQCxcuFB1gPP78vHxwYkTJ9TGXr39Oi/Py87ORmRkJCpXrqxablhYmFpICQsLg7W1NUqXLg3gRThp1KgRJk2ahLNnz8LMzAzbt2/XuC4zMzO1rUGv6+XYsWNqY2FhYahYsSJMTEzy1ZMmvr6+iIqKyjO+efNmbNu2DaGhobhz5w5++umnPHMuXboEX1/fd143UXHFkEREBbZ48WJkZ2fDz88PmzdvRnR0NGJiYrBu3TpcuXKlwGFg2LBh2Lt3L2bMmIGrV69i4cKF2Lt3b74eu2jRImzfvh1XrlzB4MGDkZycjL59+wJ48S28O3fuYOjQobhy5Qr+/vtvTJgwAaNGjYJUKsXJkycxZcoUnD59Grdv38a2bdvw8OHD135V3svLCydPnkRcXBwePXqkccvP6NGjceDAAfz000+4evUq1qxZg4ULF2LMmDEFek5e1apVK1y+fFlta9Ldu3fx9ddfY/r06WjcuDFWr16NqVOn5gmYR48eRVBQ0Hutn6hY0uUBUURkuOLj48WQIUOEt7e3kMlkokSJEqJevXpi5syZ4tmzZ6p5AMT27dvVHvvqwc5CCLFixQrh7u4uLCwsRIcOHcSsWbPydeD2hg0bxAcffCDMzMxElSpVxIEDB9TmhYaGirp16wozMzPh4uIivv32W6FQKIQQQkRFRYlWrVoJR0dHIZfLRcWKFcWCBQtUj331wO2YmBhRv359YWFhIQCI2NhYjb1s3bpV+Pj4CJlMJsqUKSNmzpypVpOnp6eYO3eu2ljNmjVV37p7nfr164ulS5cKIV4c7N2yZUvRqlUrtQPTR44cKcqVK6f6hl5YWJgoWbKkSE9Pf+OyiSgviRBaPhiAiKgIxMXFwdvbG2fPnkWtWrV0XU6R2L17N8aMGYNLly5BKs3fjoBPPvkEvr6++O677wq5OiLjwzNuExEZiLZt2+LatWu4d+8ePDw83jo/MzMTNWvWxMiRI4ugOiLjwy1JRGSQiuOWJCIqWgxJRERERBrw221EREREGjAkEREREWnAkERERESkAUMSERERkQYMSUREREQaMCQRERERacCQRERERKQBQxIRERGRBgxJRERERBr8H4SSXAZnmzXpAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# left side of the index (~max_grad_idx data)\n",
    "cut_grid = flame.grid[:max_grad_idx + 1]  # ~idx grid\n",
    "cut_u = flame.velocity[:max_grad_idx + 1]  # ~idx velocity\n",
    "\n",
    "# plot grid-velocity graph\n",
    "# plt.plot(flame.grid, flame.velocity, label=\"Velocity vs. Grid\")\n",
    "# plt.plot(cut_grid, cut_u, label=\"Velocity vs. Grid (left of max |du/dx|)\")\n",
    "plt.plot(fl1.grid, fl1.velocity, label=\"Velocity vs. Grid\")\n",
    "plt.xlabel('Grid position (x)')\n",
    "plt.ylabel('Velocity (u)')\n",
    "plt.title('Velocity Profile up to Max |du/dx|')\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "6292bdb2-00e1-428d-9733-f224c06a5da8",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1800.0\n"
     ]
    }
   ],
   "source": [
    "# Gas object 2 for the counterflow flame with non-equilibrium products\n",
    "gas2 = ct.Solution('gri30.yaml')\n",
    "gas2.TP = T0, p0\n",
    "gas2.set_equivalence_ratio(phi, fuel, oxidizer)\n",
    "\n",
    "# Create the flame simulation object\n",
    "fl2 = ct.CounterflowPremixedFlame(gas=gas2, width=width)\n",
    "fl2.set_refine_criteria(ratio=3, slope=0.1, curve=0.2, prune=0.02)\n",
    "\n",
    "# set the boundary flow rates\n",
    "fl2.reactants.mdot = mdot_reactants\n",
    "fl2.products.mdot = mdot_products\n",
    "\n",
    "#product temperature, temperature of my project\n",
    "fl2.products.T = 1800\n",
    "fl2.products.X = fl1.products.X #composition of the products\n",
    "\n",
    "fl2.set_initial_guess(equilibrate=False)\n",
    "fl2.solve(loglevel, auto=True)\n",
    "print(fl2.products.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "d89977df-308b-4861-9369-102c1cd121e4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adiabatic temperature (Tad) from FreeFlame: 1812.87 K\n"
     ]
    }
   ],
   "source": [
    "# get adiabatic temperature from the FreeFlame object\n",
    "def get_adiabatic_temperature_from_flame(flame):\n",
    "    Tad = flame.T[-1]  # final status of the flame\n",
    "    return Tad\n",
    "\n",
    "Tad = get_adiabatic_temperature_from_flame(fl1)\n",
    "print(f'Adiabatic temperature (Tad) from FreeFlame: {Tad:.2f} K')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "790733c6-d089-4538-bb81-15f612e4fd52",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "heat loss coefficient (beta): 0.992\n"
     ]
    }
   ],
   "source": [
    "def calculate_heat_loss_coefficient(Tp, Tad, Tr):\n",
    "    beta = (Tp - Tr) / (Tad - Tr)\n",
    "    return beta\n",
    "\n",
    "# set parameter value\n",
    "Tp = fl2.products.T\n",
    "Tr = T0\n",
    "\n",
    "# calculate heat loss coefficient\n",
    "beta = calculate_heat_loss_coefficient(Tp, Tad, Tr)\n",
    "print(f'heat loss coefficient (beta): {beta:.3f}')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "af6154b9-4d5f-40f3-932e-b564b3402584",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Consumption speed S_cF: 0.15381 m/s\n"
     ]
    }
   ],
   "source": [
    "def calculate_consumption_speed(flame, gas, fuel_species):\n",
    "    # unburnt density (first point of the grid)\n",
    "    rho_u = flame.density[0]  \n",
    "    \n",
    "    # calculate (Y_k,b - Y_k,u) for each fuel species\n",
    "    mass_fraction_diff = 0\n",
    "    integral_numerator = 0\n",
    "\n",
    "    # (Y_f,u, first point)\n",
    "    total_fuel_mass_fraction = 0\n",
    "    for fuel in fuel_species:\n",
    "        i_fuel = gas.species_index(fuel)\n",
    "        total_fuel_mass_fraction += flame.Y[i_fuel, 0]\n",
    "\n",
    "    for fuel in fuel_species:\n",
    "        i_fuel = gas.species_index(fuel)  # fuel index\n",
    "\n",
    "        # calculate Y_k,b (end value) and Y_k_u (first value)\n",
    "        Y_k_b = flame.Y[i_fuel, -1]  # end point (product)\n",
    "        Y_k_u = flame.Y[i_fuel, 0]   # first point (reactant)\n",
    "        mass_fraction_diff += Y_k_b - Y_k_u  # sigma (Y_k,b - Y_k,u)\n",
    "\n",
    "    # calculate nk (n_k = Y_k,u / Y_f,u)\n",
    "    ih2 = gas.species_index('H2')\n",
    "    Y_h2_u = flame.Y[ih2, 0]  # mass fraction\n",
    "    nk_h2 = Y_h2_u / total_fuel_mass_fraction\n",
    "\n",
    "    ich4 = gas.species_index('CH4')\n",
    "    Y_ch4_u = flame.Y[ich4, 0]  # mass fraction\n",
    "    nk_ch4 = Y_ch4_u / total_fuel_mass_fraction\n",
    "\n",
    "    # integral about production rates\n",
    "    integral_h2 = scipy.integrate.simpson(gas.molecular_weights[ih2] * flame.net_production_rates[ih2], x=flame.grid)\n",
    "    integral_ch4 = scipy.integrate.simpson(gas.molecular_weights[ich4] * flame.net_production_rates[ich4], x=flame.grid)\n",
    "\n",
    "    # integral_numerator update (nk*integral value)\n",
    "    integral_numerator = (nk_h2 * integral_h2) + (nk_ch4 * integral_ch4)\n",
    "\n",
    "    S_cF = integral_numerator / (rho_u * mass_fraction_diff)\n",
    "    \n",
    "    return S_cF\n",
    "\n",
    "fuel_species = ['CH4', 'H2']\n",
    "S_cF = calculate_consumption_speed(fl2, gas1, fuel_species)\n",
    "print(f'Consumption speed S_cF: {S_cF:.5f} m/s')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "d0419b13-0a57-4d4d-9854-e98c4517c11e",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[0.23629613685668938, 0.23629638445577197, 0.23629638455652693, 0.2362963841880407, 0.2362963843735158, 0.23629638416058243, 0.2362963849342203, 0.23629638446471038, 0.23629638429732028, 0.23629638424927712]\n"
     ]
    }
   ],
   "source": [
    "# Range of Ka values (logarithmic)\n",
    "Ka_values = np.logspace(-2, 2, 10)  # Logarithmic range for Ka from 10^-2 to 10^2\n",
    "S_cF_values = []\n",
    "\n",
    "for Ka in Ka_values:\n",
    "    # Calculate the new inlet velocity based on Ka\n",
    "    inlet_velocity = Ka * width / flame.density[0]  # Example equation based on Ka and width\n",
    "    flame.inlet.mdot = inlet_velocity * flame.density[0]  # 질량 유량 설정\n",
    "    flame.solve(loglevel=0, auto=True)\n",
    "\n",
    "    # Calculate consumption speed with the new boundary conditions\n",
    "    S_cF = calculate_consumption_speed(flame, gas1, fuel_species)\n",
    "    S_cF_values.append(S_cF)\n",
    "    \n",
    "# print(Ka_values)\n",
    "print(S_cF_values)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "53d22517-f4e4-4e4c-866d-773900676aaf",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plotting the results\n",
    "plt.figure()\n",
    "plt.plot(Ka_values, S_cF_values, label=r'$\\phi=0.7, X_{H_2}=0.566$', color='blue')\n",
    "plt.axhline(y=1, color='k', linestyle='--', linewidth=1)  # stagnation line\n",
    "plt.xscale('log')\n",
    "plt.title('Consumption speed for adiabatic strained laminar flame')\n",
    "plt.xlabel('Ka')\n",
    "plt.ylabel(r'$S_c/S_L^0$')\n",
    "plt.ylim(0, 2)\n",
    "plt.legend()\n",
    "plt.grid(True)\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.9"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
